#ifndef DJSS_H
#define DJSS_H

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>


#include "GadgetParticleType/gpt.01.h"
#include "Timer/timer.01.h"
#include "Vector/vector.03.h"
#include "utils/utils.02.h"

#define JERK_OFF
#define ACC_OFF
//#define POT_OFF
#include "ForceFinder/ff.01.h"

	
class PotMinFinder { 	

	public:
	
	Timer tr;
	Vector<double> PC;
	
	dj01::FJPF fjpf;

	PotMinFinder():
	potp(0)
	{}
	
	uint N;
	double *pot;
	
	
	Vector<double> find(uint ncore, uint N,  double*m, double *x0,double*x1,double*x2, double*pot,double eps2){
		
		this->N=N;
		
		fjpf.ff(N, N, m, x0,x1,x2, NULL,NULL,NULL, NULL,NULL,NULL, NULL,NULL,NULL, pot, eps2, ncore);
		
		double pot_min(1e99);
		for (int i=0; i<N; i++)
			if (pot[i]<pot_min) {
				pot_min = pot[i];
				PC.x=x0[i];
				PC.y=x1[i];
				PC.z=x2[i];
			}
		
		return PC;	
	}
	
	
	template <typename T>
	Vector<double> find(uint ncore, vector<T> &p, double eps2){
		
		tr.set();
		
		uint align = 512;
		uint N = p.size();
		uint _N = N + 2*align;
		
		char *mp, *x0p,*x1p,*x2p;
		
		double *m  = (double*)new_align(_N, 8, align, mp);
		
		double *x0 = (double*)new_align(_N, 8, align, x0p);
		double *x1 = (double*)new_align(_N, 8, align, x1p);
		double *x2 = (double*)new_align(_N, 8, align, x2p);
		
		if (potp) delete [] potp;	
		pot = (double*)new_align(_N, 8, align, potp);
		
		for (int i=0; i<N; i++)	{
			
			m[i]=p[i].m;
			x0[i]= p[i].R.x;
			x1[i]= p[i].R.y;
			x2[i]= p[i].R.z;
			
		} m[N]=0;

		PC = find(ncore, N, m, x0,x1,x2, pot,eps2);
	
		delete [] mp;
		delete [] x0p;
		delete [] x1p;
		delete [] x2p;

		tr.get();
		
		return PC;
		
	}
	
	
	
	template <typename T1,typename T2>
	Vector<T2> find(uint ncore, uint N,  T1 *_m, Vector<T2> *r, double eps2){
		
		tr.set();
		
		uint align = 512;
		uint _N = N + 2*align;
		
		char *mp, *x0p,*x1p,*x2p;
		
		double *m  = (double*)new_align(_N, 8, align, mp);
		
		double *x0 = (double*)new_align(_N, 8, align, x0p);
		double *x1 = (double*)new_align(_N, 8, align, x1p);
		double *x2 = (double*)new_align(_N, 8, align, x2p);
		
		if (potp) delete [] potp;	
		pot = (double*)new_align(_N, 8, align, potp);
		
		for (int i=0; i<N; i++)	{
			
			m[i]=_m[i];
			x0[i]= r[i].x;
			x1[i]= r[i].y;
			x2[i]= r[i].z;
			
		} m[N]=0;

		PC = find(ncore, N, m, x0,x1,x2, pot,eps2);
	
		delete [] mp;
		delete [] x0p;
		delete [] x1p;
		delete [] x2p;

		tr.get();
		
		return PC;
		
	}
	
	
	template <typename T1,typename T2, typename T3>
	Vector<T2> find(uint ncore, uint N,  T1 *p, T2 *_m, T3 *_r, double eps2){

		tr.set();
		
		uint align = 512;
		uint _N = N + 2*align;
		
		char *mp, *x0p,*x1p,*x2p;
		
		double *m  = (double*)new_align(_N, 8, align, mp);
		double *x0 = (double*)new_align(_N, 8, align, x0p);
		double *x1 = (double*)new_align(_N, 8, align, x1p);
		double *x2 = (double*)new_align(_N, 8, align, x2p);
		
		if (potp) delete [] potp;	
		pot = (double*)new_align(_N, 8, align, potp);
		
		for (int i=0; i<N; i++)	{
			
			m[i]= *(T2*)&((T1*)_m)[i];
			x0[i]= (*(T3*)&((T1*)_r)[i]).x;
			x1[i]= (*(T3*)&((T1*)_r)[i]).y;
			x2[i]= (*(T3*)&((T1*)_r)[i]).z;
			
		} m[N]=0;

		PC = find(ncore, N, m, x0,x1,x2, pot,eps2);
	
		delete [] mp;
		delete [] x0p;
		delete [] x1p;
		delete [] x2p;

		tr.get();
		
		return PC;
		
		
		return 0;
	}

	// in Vector in vector container
	template <typename T1,typename T2>
	Vector<double> find(uint ncore, const vector<T1> &_m, const vector< Vector<T2> > &r, double eps2){
		
		tr.set();
		
		uint align = 512;
		N = r.size();
		uint _N = N + 2*align;
		
		char *mp, *x0p,*x1p,*x2p;
		
		double *m  = (double*)new_align(_N, 8, align, mp);
		
		double *x0 = (double*)new_align(_N, 8, align, x0p);
		double *x1 = (double*)new_align(_N, 8, align, x1p);
		double *x2 = (double*)new_align(_N, 8, align, x2p);
		
		if (potp) delete [] potp;	
		pot = (double*)new_align(_N, 8, align, potp);
		
		for (int i=0; i<N; i++)	{
			
			m[i]=_m[i];
			x0[i]= r[i].x;
			x1[i]= r[i].y;
			x2[i]= r[i].z;
			
		} m[N]=0;

		PC = find(ncore, N, m, x0,x1,x2, pot,eps2);
	
		delete [] mp;
		delete [] x0p;
		delete [] x1p;
		delete [] x2p;

		tr.get();
		
		return PC;
		
	}


	char *potp;
	
	
	Vector<double> find(uint ncore, uint ofs, uint N,  P *p, double eps2){
		
		tr.set();
		
		uint align = 512;
		uint _N = N + 2*align;
		
		char *mp, *x0p,*x1p,*x2p;
		
		double *m  = (double*)new_align(_N, 8, align, mp);
		
		double *x0 = (double*)new_align(_N, 8, align, x0p);
		double *x1 = (double*)new_align(_N, 8, align, x1p);
		double *x2 = (double*)new_align(_N, 8, align, x2p);
		
		if (potp) delete [] potp;
		pot = (double*)new_align(_N, 8, align, potp);
		
		for (int i=0; i<N; i++)	{
			
			m[i] = p[ofs+i].Mass;
			x0[i]= p[ofs+i].Pos[0];
			x1[i]= p[ofs+i].Pos[1];
			x2[i]= p[ofs+i].Pos[2];
			
		} m[N]=0;
		

		PC = find(ncore, N, m, x0,x1,x2, pot,eps2);
		
		delete [] mp;
		delete [] x0p;
		delete [] x1p;
		delete [] x2p;
		
				
		tr.get();
		
		return PC;

		
	}
			
	void print(){

		printf("\033[22;34mProcessingTime:\033[22;30m %g\n", tr.sin);
		printf("\033[22;34mPotential Center:\033[22;30m %g, %g, %g\n", PC.x, PC.y, PC.z);
		
	}
	
};

typedef PotMinFinder PMF;

#endif
